Clean data

Load data

dat <- readxl::read_xlsx("/Volumes/manhnd/Yi_roe_data.xlsx")
# filter variables for analyis
dat <- dat %>% dplyr::select("serialno", "country", "papm", "gender", "age", "education", "marital", "occupation", "live_area", "mobile", "internet")

Transform data

# Remove the gender as others (18/2912)
dat1 <- dat %>%
  filter(gender != 3 & gender != 4) %>%
  mutate(marital = dplyr::recode(marital, `4` = 3, `5` = 3))

# Transform the valuables
dat1 <- dat1 %>%
  mutate(
    papm = factor(papm, ordered = TRUE),
    gender = factor(gender, labels = c("Male", "Female")),
    marital = factor(marital, labels = c("Never married", "Married/living with a partner", "Separated/widowed/divorced")),
    education = factor(education, labels = c("No formal education", "Primary school", "Secondary school", "Diploma", "Degree")),
    occupation = factor(occupation, labels = c("Employed", "Unemployed", "Student", "Retiree", "Homemaker")),
    live_area = factor(live_area, labels = c("Urban", "Peri-urban", "Rural", "Slums")),
    mobile = factor(mobile, labels = c("No", "Yes feature phone", "Yes smart phone")),
    internet = factor(internet, labels = c("No", "Yes")),
    country = factor(country)
  )

summary(dat1)
##     serialno             country    papm        gender          age       
##  Min.   :   1.0   Philippines:459   1:1596   Male  :1365   Min.   :18.00  
##  1st Qu.: 726.2   Uganda     :395   2: 291   Female:1529   1st Qu.:22.00  
##  Median :1453.5   Indonesia  :368   3: 285                 Median :30.00  
##  Mean   :1455.6   Zimbabwe   :292   4: 179                 Mean   :34.95  
##  3rd Qu.:2185.8   Cameroon   :289   5: 543                 3rd Qu.:44.00  
##  Max.   :2912.0   Bangladesh :280                          Max.   :90.00  
##                   (Other)    :811                                         
##                education                             marital    
##  No formal education: 173   Never married                :1221  
##  Primary school     : 476   Married/living with a partner:1395  
##  Secondary school   :1182   Separated/widowed/divorced   : 278  
##  Diploma            : 405                                       
##  Degree             : 658                                       
##                                                                 
##                                                                 
##       occupation        live_area                  mobile     internet  
##  Employed  :1294   Urban     :1066   No               : 257   No : 839  
##  Unemployed: 566   Peri-urban: 374   Yes feature phone: 626   Yes:2055  
##  Student   : 681   Rural     :1281   Yes smart phone  :2011             
##  Retiree   :  65   Slums     : 173                                      
##  Homemaker : 288                                                        
##                                                                         
## 
str(dat1)
## tibble [2,894 × 11] (S3: tbl_df/tbl/data.frame)
##  $ serialno  : num [1:2894] 1 2 3 4 5 6 7 8 9 10 ...
##  $ country   : Factor w/ 9 levels "Bangladesh","Cameroon",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ papm      : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ gender    : Factor w/ 2 levels "Male","Female": 2 2 2 2 2 1 1 2 1 2 ...
##  $ age       : num [1:2894] 29 36 52 27 31 31 70 50 49 45 ...
##  $ education : Factor w/ 5 levels "No formal education",..: 2 1 1 2 2 2 1 1 1 1 ...
##  $ marital   : Factor w/ 3 levels "Never married",..: 2 2 3 2 2 2 2 2 2 2 ...
##  $ occupation: Factor w/ 5 levels "Employed","Unemployed",..: 5 5 1 5 5 1 2 5 1 5 ...
##  $ live_area : Factor w/ 4 levels "Urban","Peri-urban",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ mobile    : Factor w/ 3 levels "No","Yes feature phone",..: 2 2 2 2 2 2 2 1 2 2 ...
##  $ internet  : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...

Visualize some variables

papm versus countries

dat2 <- dat1 %>%
  select("country", "papm") %>%
  group_by(country, papm) %>%
  summarise(n = n()) %>%
  mutate(proportion_papm = n / sum(n)) %>%
  ungroup()
## `summarise()` has grouped output by 'country'. You can override using the
## `.groups` argument.
ggplot(data = dat2, aes(x = country, y = proportion_papm)) +
  geom_col(width = 0.5) +
  facet_grid(~papm) +
  coord_flip()

Education versus countries

dat2 <- dat %>%
  select("country", "education") %>%
  group_by(country, education) %>%
  summarise(n = n()) %>%
  mutate(proportion_education = n / sum(n)) %>%
  ungroup()
## `summarise()` has grouped output by 'country'. You can override using the
## `.groups` argument.
ggplot(data = dat2, aes(x = country, y = proportion_education)) +
  geom_col(width = 0.5) +
  facet_grid(~education) +
  coord_flip()

MASS package: polr modell

Fit the modell without country interaction

m1 <- polr(formula = papm ~ country + gender + age + education + marital + occupation + live_area + mobile + internet, data = dat1, Hess = TRUE)

summary(m1)
## Call:
## polr(formula = papm ~ country + gender + age + education + marital + 
##     occupation + live_area + mobile + internet, data = dat1, 
##     Hess = TRUE)
## 
## Coefficients:
##                                          Value Std. Error  t value
## countryCameroon                       1.030800   0.304414  3.38618
## countryIndia                          3.224494   0.272151 11.84819
## countryIndonesia                      1.723045   0.279535  6.16397
## countryKenya                          2.155655   0.310191  6.94944
## countryNepal                          3.933541   0.282818 13.90837
## countryPhilippines                    1.061534   0.286374  3.70681
## countryUganda                         4.131868   0.272311 15.17337
## countryZimbabwe                       4.363371   0.287158 15.19501
## genderFemale                         -0.005772   0.085817 -0.06726
## age                                  -0.001336   0.003866 -0.34557
## educationPrimary school               0.480407   0.190301  2.52446
## educationSecondary school             0.858222   0.195790  4.38338
## educationDiploma                      1.083491   0.222924  4.86036
## educationDegree                       1.204571   0.220266  5.46872
## maritalMarried/living with a partner  0.243062   0.122772  1.97978
## maritalSeparated/widowed/divorced     0.602561   0.181871  3.31312
## occupationUnemployed                 -0.183360   0.118656 -1.54530
## occupationStudent                    -0.156333   0.134595 -1.16151
## occupationRetiree                    -0.556455   0.284433 -1.95637
## occupationHomemaker                  -0.404096   0.151572 -2.66603
## live_areaPeri-urban                  -0.019113   0.135460 -0.14110
## live_areaRural                        0.056872   0.116823  0.48682
## live_areaSlums                        0.341699   0.225495  1.51533
## mobileYes feature phone               0.576591   0.161966  3.55994
## mobileYes smart phone                 0.217389   0.212906  1.02106
## internetYes                           0.491155   0.168225  2.91964
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2  4.2778  0.3845    11.1264
## 2|3  4.9192  0.3873    12.7018
## 3|4  5.5729  0.3902    14.2807
## 4|5  6.0354  0.3921    15.3938
## 
## Residual Deviance: 6159.317 
## AIC: 6219.317

Fit the model with country interaction

# First fit the model with interaction between country and all of the independent variables but the model cannot run. Then I removed the interaction between live_area and country, model can run as below
m2 <- polr(formula = papm ~ country + gender + age + education + marital + occupation + live_area + mobile + internet + gender:country + age:country + marital:country + occupation:country + mobile:country + internet:country + education:country, data = dat1, Hess = TRUE)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(m2)
## Call:
## polr(formula = papm ~ country + gender + age + education + marital + 
##     occupation + live_area + mobile + internet + gender:country + 
##     age:country + marital:country + occupation:country + mobile:country + 
##     internet:country + education:country, data = dat1, Hess = TRUE)
## 
## Coefficients:
##                                                              Value Std. Error
## countryCameroon                                          14.609312  1.400e+00
## countryIndia                                             24.913158  1.278e+00
## countryIndonesia                                         -3.055593  1.099e+00
## countryKenya                                             26.698692  1.697e+00
## countryNepal                                             27.589432  1.359e+00
## countryPhilippines                                       23.911632  1.497e+00
## countryUganda                                            26.874774  1.138e+00
## countryZimbabwe                                          16.276852  1.281e+00
## genderFemale                                              0.354487  5.578e-01
## age                                                       0.077100  3.809e-02
## educationPrimary school                                  -4.080201  3.335e-01
## educationSecondary school                                 9.866479  7.889e-01
## educationDiploma                                         11.217361  6.772e-01
## educationDegree                                          11.537265  5.694e-01
## maritalMarried/living with a partner                     -0.857760  1.184e+00
## maritalSeparated/widowed/divorced                        -2.453037  2.089e-01
## occupationUnemployed                                    -16.542871  1.461e-01
## occupationStudent                                         0.048204  1.308e+00
## occupationRetiree                                       -18.167315  3.829e-01
## occupationHomemaker                                     -12.483776  2.274e-01
## live_areaPeri-urban                                      -0.090668  1.443e-01
## live_areaRural                                           -0.001955  1.294e-01
## live_areaSlums                                            0.249086  2.547e-01
## mobileYes feature phone                                  11.977607  8.830e-01
## mobileYes smart phone                                    11.867155  1.091e+00
## internetYes                                               0.097749  1.624e+00
## countryCameroon:genderFemale                             -0.475020  6.699e-01
## countryIndia:genderFemale                                 0.074489  6.111e-01
## countryIndonesia:genderFemale                            -0.611408  6.207e-01
## countryKenya:genderFemale                                -0.179755  6.525e-01
## countryNepal:genderFemale                                -0.674883  6.148e-01
## countryPhilippines:genderFemale                          -0.091360  6.396e-01
## countryUganda:genderFemale                               -0.477477  5.858e-01
## countryZimbabwe:genderFemale                             -0.918931  6.375e-01
## countryCameroon:age                                      -0.074785  4.493e-02
## countryIndia:age                                         -0.063382  3.976e-02
## countryIndonesia:age                                     -0.033517  4.247e-02
## countryKenya:age                                         -0.031155  4.203e-02
## countryNepal:age                                         -0.064597  4.016e-02
## countryPhilippines:age                                   -0.080832  4.119e-02
## countryUganda:age                                        -0.073202  3.865e-02
## countryZimbabwe:age                                      -0.122026  4.188e-02
## countryCameroon:maritalMarried/living with a partner      1.420520  1.273e+00
## countryIndia:maritalMarried/living with a partner         0.976123  1.236e+00
## countryIndonesia:maritalMarried/living with a partner     0.718878  1.232e+00
## countryKenya:maritalMarried/living with a partner         0.783927  1.265e+00
## countryNepal:maritalMarried/living with a partner         0.381139  1.243e+00
## countryPhilippines:maritalMarried/living with a partner   1.822441  1.282e+00
## countryUganda:maritalMarried/living with a partner        0.919528  1.221e+00
## countryZimbabwe:maritalMarried/living with a partner      2.381148  1.259e+00
## countryCameroon:maritalSeparated/widowed/divorced       -12.787097  8.655e-08
## countryIndia:maritalSeparated/widowed/divorced            2.628711  6.726e-01
## countryIndonesia:maritalSeparated/widowed/divorced        0.975907  6.736e-01
## countryKenya:maritalSeparated/widowed/divorced            2.709081  5.056e-01
## countryNepal:maritalSeparated/widowed/divorced            3.123797  6.262e-01
## countryPhilippines:maritalSeparated/widowed/divorced      3.826500  6.370e-01
## countryUganda:maritalSeparated/widowed/divorced           2.576574  4.045e-01
## countryZimbabwe:maritalSeparated/widowed/divorced         4.439632  5.209e-01
## countryCameroon:occupationUnemployed                     16.766464  4.480e-01
## countryIndia:occupationUnemployed                        16.584719  4.851e-01
## countryIndonesia:occupationUnemployed                    15.559812  4.244e-01
## countryKenya:occupationUnemployed                        16.774798  4.632e-01
## countryNepal:occupationUnemployed                        14.736578  5.290e-01
## countryPhilippines:occupationUnemployed                  16.481324  4.856e-01
## countryUganda:occupationUnemployed                       16.816579  2.401e-01
## countryZimbabwe:occupationUnemployed                     15.691493  3.303e-01
## countryCameroon:occupationStudent                        -0.162487  1.393e+00
## countryIndia:occupationStudent                            0.259397  1.358e+00
## countryIndonesia:occupationStudent                       -0.397042  1.361e+00
## countryKenya:occupationStudent                           -0.375236  1.443e+00
## countryNepal:occupationStudent                           -0.461873  1.359e+00
## countryPhilippines:occupationStudent                      0.709372  1.412e+00
## countryUganda:occupationStudent                           0.011918  1.360e+00
## countryZimbabwe:occupationStudent                         0.052282  1.394e+00
## countryCameroon:occupationRetiree                         0.304758  1.344e-08
## countryIndia:occupationRetiree                           17.765439  6.078e-01
## countryIndonesia:occupationRetiree                        1.689434  1.391e-07
## countryKenya:occupationRetiree                           18.059350  1.465e+00
## countryNepal:occupationRetiree                           16.820388  5.936e-01
## countryPhilippines:occupationRetiree                     18.007567  8.731e-01
## countryUganda:occupationRetiree                          17.472933  7.974e-01
## countryZimbabwe:occupationRetiree                        36.326893  6.540e-08
## countryCameroon:occupationHomemaker                      11.750979  9.985e-01
## countryIndia:occupationHomemaker                         11.743691  3.820e-01
## countryIndonesia:occupationHomemaker                     11.075414  6.715e-01
## countryKenya:occupationHomemaker                         13.993172  7.158e-01
## countryNepal:occupationHomemaker                         11.784830  4.632e-01
## countryPhilippines:occupationHomemaker                   12.831394  6.032e-01
## countryUganda:occupationHomemaker                        12.899331  3.283e-01
## countryZimbabwe:occupationHomemaker                      10.611018  9.371e-01
## countryCameroon:mobileYes feature phone                 -12.848847  1.267e+00
## countryIndia:mobileYes feature phone                    -11.411442  9.488e-01
## countryIndonesia:mobileYes feature phone                -14.822680  2.206e+00
## countryKenya:mobileYes feature phone                    -11.090439  1.084e+00
## countryNepal:mobileYes feature phone                    -12.424780  1.400e+00
## countryPhilippines:mobileYes feature phone              -11.984637  1.099e+00
## countryUganda:mobileYes feature phone                   -11.399280  9.077e-01
## countryZimbabwe:mobileYes feature phone                  -8.440743  1.168e+00
## countryCameroon:mobileYes smart phone                   -11.923094  1.571e+00
## countryIndia:mobileYes smart phone                      -11.281882  1.357e+00
## countryIndonesia:mobileYes smart phone                  -13.735593  2.200e+00
## countryKenya:mobileYes smart phone                      -10.639113  1.271e+00
## countryNepal:mobileYes smart phone                      -12.697334  1.963e+00
## countryPhilippines:mobileYes smart phone                -13.246191  1.279e+00
## countryUganda:mobileYes smart phone                     -11.300390  1.198e+00
## countryZimbabwe:mobileYes smart phone                   -10.451132  1.268e+00
## countryCameroon:internetYes                              -0.142093  1.898e+00
## countryIndia:internetYes                                  0.343617  1.800e+00
## countryIndonesia:internetYes                             19.075766  1.099e+00
## countryKenya:internetYes                                 -0.298728  1.683e+00
## countryNepal:internetYes                                  1.112136  2.262e+00
## countryPhilippines:internetYes                            2.030549  1.746e+00
## countryUganda:internetYes                                 0.590445  1.694e+00
## countryZimbabwe:internetYes                              -0.705238  1.672e+00
## countryCameroon:educationPrimary school                  13.077706  6.222e-01
## countryIndia:educationPrimary school                      4.762476  5.766e-01
## countryIndonesia:educationPrimary school                 16.204295  7.376e-01
## countryKenya:educationPrimary school                      0.185179  1.174e+00
## countryNepal:educationPrimary school                      3.920423  5.980e-01
## countryPhilippines:educationPrimary school                3.667566  8.605e-01
## countryUganda:educationPrimary school                     4.581924  3.971e-01
## countryZimbabwe:educationPrimary school                  17.735081  7.247e-01
## countryCameroon:educationSecondary school                 0.958803  8.201e-01
## countryIndia:educationSecondary school                   -8.921375  9.204e-01
## countryIndonesia:educationSecondary school                1.257021  7.397e-01
## countryKenya:educationSecondary school                  -12.508199  1.331e+00
## countryNepal:educationSecondary school                   -9.343525  9.348e-01
## countryPhilippines:educationSecondary school             -9.735912  1.081e+00
## countryUganda:educationSecondary school                  -9.409909  8.247e-01
## countryZimbabwe:educationSecondary school                 3.027823  7.846e-01
## countryCameroon:educationDiploma                          0.244220  7.472e-01
## countryIndia:educationDiploma                            -9.866073  8.812e-01
## countryIndonesia:educationDiploma                         0.688810  6.993e-01
## countryKenya:educationDiploma                           -14.272035  1.290e+00
## countryNepal:educationDiploma                            -9.862474  8.626e-01
## countryPhilippines:educationDiploma                     -11.230641  1.041e+00
## countryUganda:educationDiploma                          -10.547329  7.987e-01
## countryZimbabwe:educationDiploma                          1.256645  7.247e-01
## countryCameroon:educationDegree                          -0.181528  6.094e-01
## countryIndia:educationDegree                             -9.934391  7.724e-01
## countryIndonesia:educationDegree                          0.151529  5.295e-01
## countryKenya:educationDegree                            -13.209072  1.288e+00
## countryNepal:educationDegree                            -10.477752  7.820e-01
## countryPhilippines:educationDegree                      -11.690755  9.826e-01
## countryUganda:educationDegree                           -11.004074  1.124e+00
## countryZimbabwe:educationDegree                           0.496835  5.932e-01
##                                                            t value
## countryCameroon                                          1.043e+01
## countryIndia                                             1.950e+01
## countryIndonesia                                        -2.779e+00
## countryKenya                                             1.574e+01
## countryNepal                                             2.030e+01
## countryPhilippines                                       1.597e+01
## countryUganda                                            2.362e+01
## countryZimbabwe                                          1.270e+01
## genderFemale                                             6.355e-01
## age                                                      2.024e+00
## educationPrimary school                                 -1.224e+01
## educationSecondary school                                1.251e+01
## educationDiploma                                         1.656e+01
## educationDegree                                          2.026e+01
## maritalMarried/living with a partner                    -7.245e-01
## maritalSeparated/widowed/divorced                       -1.174e+01
## occupationUnemployed                                    -1.132e+02
## occupationStudent                                        3.684e-02
## occupationRetiree                                       -4.744e+01
## occupationHomemaker                                     -5.490e+01
## live_areaPeri-urban                                     -6.285e-01
## live_areaRural                                          -1.511e-02
## live_areaSlums                                           9.779e-01
## mobileYes feature phone                                  1.356e+01
## mobileYes smart phone                                    1.088e+01
## internetYes                                              6.018e-02
## countryCameroon:genderFemale                            -7.091e-01
## countryIndia:genderFemale                                1.219e-01
## countryIndonesia:genderFemale                           -9.851e-01
## countryKenya:genderFemale                               -2.755e-01
## countryNepal:genderFemale                               -1.098e+00
## countryPhilippines:genderFemale                         -1.428e-01
## countryUganda:genderFemale                              -8.151e-01
## countryZimbabwe:genderFemale                            -1.442e+00
## countryCameroon:age                                     -1.664e+00
## countryIndia:age                                        -1.594e+00
## countryIndonesia:age                                    -7.893e-01
## countryKenya:age                                        -7.413e-01
## countryNepal:age                                        -1.609e+00
## countryPhilippines:age                                  -1.962e+00
## countryUganda:age                                       -1.894e+00
## countryZimbabwe:age                                     -2.914e+00
## countryCameroon:maritalMarried/living with a partner     1.116e+00
## countryIndia:maritalMarried/living with a partner        7.898e-01
## countryIndonesia:maritalMarried/living with a partner    5.834e-01
## countryKenya:maritalMarried/living with a partner        6.198e-01
## countryNepal:maritalMarried/living with a partner        3.066e-01
## countryPhilippines:maritalMarried/living with a partner  1.421e+00
## countryUganda:maritalMarried/living with a partner       7.534e-01
## countryZimbabwe:maritalMarried/living with a partner     1.892e+00
## countryCameroon:maritalSeparated/widowed/divorced       -1.477e+08
## countryIndia:maritalSeparated/widowed/divorced           3.908e+00
## countryIndonesia:maritalSeparated/widowed/divorced       1.449e+00
## countryKenya:maritalSeparated/widowed/divorced           5.359e+00
## countryNepal:maritalSeparated/widowed/divorced           4.989e+00
## countryPhilippines:maritalSeparated/widowed/divorced     6.007e+00
## countryUganda:maritalSeparated/widowed/divorced          6.369e+00
## countryZimbabwe:maritalSeparated/widowed/divorced        8.523e+00
## countryCameroon:occupationUnemployed                     3.743e+01
## countryIndia:occupationUnemployed                        3.419e+01
## countryIndonesia:occupationUnemployed                    3.666e+01
## countryKenya:occupationUnemployed                        3.622e+01
## countryNepal:occupationUnemployed                        2.786e+01
## countryPhilippines:occupationUnemployed                  3.394e+01
## countryUganda:occupationUnemployed                       7.003e+01
## countryZimbabwe:occupationUnemployed                     4.750e+01
## countryCameroon:occupationStudent                       -1.166e-01
## countryIndia:occupationStudent                           1.910e-01
## countryIndonesia:occupationStudent                      -2.918e-01
## countryKenya:occupationStudent                          -2.601e-01
## countryNepal:occupationStudent                          -3.400e-01
## countryPhilippines:occupationStudent                     5.024e-01
## countryUganda:occupationStudent                          8.762e-03
## countryZimbabwe:occupationStudent                        3.751e-02
## countryCameroon:occupationRetiree                        2.267e+07
## countryIndia:occupationRetiree                           2.923e+01
## countryIndonesia:occupationRetiree                       1.215e+07
## countryKenya:occupationRetiree                           1.233e+01
## countryNepal:occupationRetiree                           2.834e+01
## countryPhilippines:occupationRetiree                     2.063e+01
## countryUganda:occupationRetiree                          2.191e+01
## countryZimbabwe:occupationRetiree                        5.555e+08
## countryCameroon:occupationHomemaker                      1.177e+01
## countryIndia:occupationHomemaker                         3.075e+01
## countryIndonesia:occupationHomemaker                     1.649e+01
## countryKenya:occupationHomemaker                         1.955e+01
## countryNepal:occupationHomemaker                         2.544e+01
## countryPhilippines:occupationHomemaker                   2.127e+01
## countryUganda:occupationHomemaker                        3.929e+01
## countryZimbabwe:occupationHomemaker                      1.132e+01
## countryCameroon:mobileYes feature phone                 -1.014e+01
## countryIndia:mobileYes feature phone                    -1.203e+01
## countryIndonesia:mobileYes feature phone                -6.720e+00
## countryKenya:mobileYes feature phone                    -1.023e+01
## countryNepal:mobileYes feature phone                    -8.874e+00
## countryPhilippines:mobileYes feature phone              -1.090e+01
## countryUganda:mobileYes feature phone                   -1.256e+01
## countryZimbabwe:mobileYes feature phone                 -7.224e+00
## countryCameroon:mobileYes smart phone                   -7.591e+00
## countryIndia:mobileYes smart phone                      -8.314e+00
## countryIndonesia:mobileYes smart phone                  -6.244e+00
## countryKenya:mobileYes smart phone                      -8.372e+00
## countryNepal:mobileYes smart phone                      -6.467e+00
## countryPhilippines:mobileYes smart phone                -1.035e+01
## countryUganda:mobileYes smart phone                     -9.429e+00
## countryZimbabwe:mobileYes smart phone                   -8.240e+00
## countryCameroon:internetYes                             -7.487e-02
## countryIndia:internetYes                                 1.909e-01
## countryIndonesia:internetYes                             1.735e+01
## countryKenya:internetYes                                -1.775e-01
## countryNepal:internetYes                                 4.916e-01
## countryPhilippines:internetYes                           1.163e+00
## countryUganda:internetYes                                3.485e-01
## countryZimbabwe:internetYes                             -4.219e-01
## countryCameroon:educationPrimary school                  2.102e+01
## countryIndia:educationPrimary school                     8.260e+00
## countryIndonesia:educationPrimary school                 2.197e+01
## countryKenya:educationPrimary school                     1.578e-01
## countryNepal:educationPrimary school                     6.556e+00
## countryPhilippines:educationPrimary school               4.262e+00
## countryUganda:educationPrimary school                    1.154e+01
## countryZimbabwe:educationPrimary school                  2.447e+01
## countryCameroon:educationSecondary school                1.169e+00
## countryIndia:educationSecondary school                  -9.693e+00
## countryIndonesia:educationSecondary school               1.699e+00
## countryKenya:educationSecondary school                  -9.395e+00
## countryNepal:educationSecondary school                  -9.995e+00
## countryPhilippines:educationSecondary school            -9.008e+00
## countryUganda:educationSecondary school                 -1.141e+01
## countryZimbabwe:educationSecondary school                3.859e+00
## countryCameroon:educationDiploma                         3.269e-01
## countryIndia:educationDiploma                           -1.120e+01
## countryIndonesia:educationDiploma                        9.850e-01
## countryKenya:educationDiploma                           -1.107e+01
## countryNepal:educationDiploma                           -1.143e+01
## countryPhilippines:educationDiploma                     -1.079e+01
## countryUganda:educationDiploma                          -1.321e+01
## countryZimbabwe:educationDiploma                         1.734e+00
## countryCameroon:educationDegree                         -2.979e-01
## countryIndia:educationDegree                            -1.286e+01
## countryIndonesia:educationDegree                         2.862e-01
## countryKenya:educationDegree                            -1.025e+01
## countryNepal:educationDegree                            -1.340e+01
## countryPhilippines:educationDegree                      -1.190e+01
## countryUganda:educationDegree                           -9.794e+00
## countryZimbabwe:educationDegree                          8.376e-01
## 
## Intercepts:
##     Value         Std. Error    t value      
## 1|2  2.706380e+01  1.076900e+00  2.513160e+01
## 2|3  2.775100e+01  1.077600e+00  2.575330e+01
## 3|4  2.845470e+01  1.078200e+00  2.638990e+01
## 4|5  2.895350e+01  1.078700e+00  2.684100e+01
## 
## Residual Deviance: 5866.142 
## AIC: 6166.142
# compare with model 1
anova(m1, m2, test = "Chisq")
## Likelihood ratio tests of ordinal regression models
## 
## Response: papm
##                                                                                                                                                                                                                     Model
## 1                                                                                                                               country + gender + age + education + marital + occupation + live_area + mobile + internet
## 2 country + gender + age + education + marital + occupation + live_area + mobile + internet + gender:country + age:country + marital:country + occupation:country + mobile:country + internet:country + education:country
##   Resid. df Resid. Dev   Test    Df LR stat.      Pr(Chi)
## 1      2864   6159.317                                   
## 2      2744   5866.142 1 vs 2   120 293.1756 1.110223e-16
Anova(m2)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Analysis of Deviance Table (Type II tests)
## 
## Response: papm
##                    LR Chisq Df Pr(>Chisq)    
## country             1052.99  8  < 2.2e-16 ***
## gender                 0.41  1  0.5230524    
## age                    3.39  1  0.0654129 .  
## education             28.66  4  9.161e-06 ***
## marital                3.83  2  0.1470023    
## occupation             7.84  4  0.0975951 .  
## live_area              1.64  3  0.6506180    
## mobile                15.12  2  0.0005211 ***
## internet               2.95  1  0.0860821 .  
## country:gender        10.36  8  0.2405672    
## country:age           22.71  8  0.0037536 ** 
## country:marital       31.38 16  0.0120169 *  
## country:occupation    59.74 32  0.0020828 ** 
## country:mobile        41.60 16  0.0004521 ***
## country:internet      28.18  8  0.0004409 ***
## country:education     70.05 32  0.0001168 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# m2 is beter than m1.

Testing the proportional odds assumption

https://peopleanalytics-regression-book.org/ord-reg.html:

dat_test1 <- dat1 %>% mutate(
  test1 = ifelse(papm == "1", 0, 1),
  test2 = ifelse(papm == "1" | papm == "2", 0, 1),
  test3 = ifelse(papm == "1" | papm == "2" | papm == "3", 0, 1),
  test4 = ifelse(papm == "1" | papm == "2" | papm == "3" | papm == "4", 0, 1)
)
table(dat_test1$papm)
## 
##    1    2    3    4    5 
## 1596  291  285  179  543
m1 <- glm(test1 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
  data = dat_test1 %>% filter(country == "Bangladesh"),
  family = "binomial"
)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
m2 <- glm(test2 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
  data = dat_test1 %>% filter(country == "Bangladesh"),
  family = "binomial"
)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
m3 <- glm(test3 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
  data = dat_test1 %>% filter(country == "Bangladesh"),
  family = "binomial"
)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
m4 <- glm(test4 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
  data = dat_test1 %>% filter(country == "Bangladesh"),
  family = "binomial"
)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
coefficient_comparison <- data.frame(
  test1 = summary(m1)$coefficients[, "Estimate"],
  test2 = summary(m2)$coefficients[, "Estimate"],
  test3 = summary(m3)$coefficients[, "Estimate"],
  test4 = summary(m4)$coefficients[, "Estimate"],
  diff1_2 = summary(m2)$coefficients[, "Estimate"] - summary(m1)$coefficients[, "Estimate"],
  diff1_3 = summary(m3)$coefficients[, "Estimate"] - summary(m1)$coefficients[, "Estimate"],
  diff1_4 = summary(m4)$coefficients[, "Estimate"] - summary(m1)$coefficients[, "Estimate"]
)
coefficient_comparison
##                                             test1        test2        test3
## (Intercept)                          -38.43425101 -38.43425101 -40.93666714
## genderFemale                           0.38700242   0.38700242   0.09338424
## age                                    0.08512152   0.08512152   0.08291879
## educationPrimary school                1.53180401   1.53180401   1.49493880
## educationSecondary school             17.77258036  17.77258036  18.95926562
## educationDiploma                      19.23192113  19.23192113  19.17753807
## educationDegree                       19.54783604  19.54783604  20.62539567
## maritalMarried/living with a partner  -0.65031291  -0.65031291  -0.56711481
## maritalSeparated/widowed/divorced     -2.11973287  -2.11973287  -1.11113159
## occupationUnemployed                 -18.58084035 -18.58084035 -19.58069143
## occupationStudent                      0.11277406   0.11277406  -0.67586541
## occupationRetiree                    -20.44333003 -20.44333003 -20.16806042
## occupationHomemaker                  -16.21747376 -16.21747376 -16.88829039
## live_areaPeri-urban                   -0.45811434  -0.45811434   1.00093169
## live_areaRural                         0.77632174   0.77632174   2.28880765
## mobileYes feature phone               14.16567476  14.16567476  14.22825219
## mobileYes smart phone                 14.11045566  14.11045566  14.20706684
## internetYes                            0.28136377   0.28136377   0.38110277
##                                             test4 diff1_2      diff1_3
## (Intercept)                          -61.18243089       0 -2.502416133
## genderFemale                          -0.11672626       0 -0.293618176
## age                                    0.03502682       0 -0.002202736
## educationPrimary school                1.07945192       0 -0.036865204
## educationSecondary school             19.16571382       0  1.186685264
## educationDiploma                      -2.44021033       0 -0.054383058
## educationDegree                       20.91176452       0  1.077559635
## maritalMarried/living with a partner  17.96932349       0  0.083198108
## maritalSeparated/widowed/divorced     18.74805329       0  1.008601275
## occupationUnemployed                 -16.15798504       0 -0.999851084
## occupationStudent                     -9.07947685       0 -0.788639469
## occupationRetiree                     -1.73019043       0  0.275269612
## occupationHomemaker                  -17.96936153       0 -0.670816628
## live_areaPeri-urban                   18.80811464       0  1.459046033
## live_areaRural                        20.42660607       0  1.512485918
## mobileYes feature phone               -0.59054383       0  0.062577434
## mobileYes smart phone                  5.33094271       0  0.096611175
## internetYes                           -5.93206007       0  0.099739006
##                                           diff1_4
## (Intercept)                          -22.74817988
## genderFemale                          -0.50372868
## age                                   -0.05009471
## educationPrimary school               -0.45235209
## educationSecondary school              1.39313347
## educationDiploma                     -21.67213146
## educationDegree                        1.36392848
## maritalMarried/living with a partner  18.61963640
## maritalSeparated/widowed/divorced     20.86778615
## occupationUnemployed                   2.42285531
## occupationStudent                     -9.19225091
## occupationRetiree                     18.71313961
## occupationHomemaker                   -1.75188777
## live_areaPeri-urban                   19.26622898
## live_areaRural                        19.65028434
## mobileYes feature phone              -14.75621859
## mobileYes smart phone                 -8.77951296
## internetYes                           -6.21342384

We could see there is a large difference coefficients. The assumption of proportional odds is violated. We need to use different methods. The book recommend some other methods explain more details in Agresti (2010))

  • Baseline logistic model. This model is the same as the multinomial regression model covered in the previous chapter, using the lowest ordinal value as the reference.
  • Adjacent-category logistic model. This model compares each level of the ordinal variable to the next highest level, and it is a constrained version of the baseline logistic model. The brglm2 package in R offers a function bracl() for calculating an adjacent category logistic model.
  • Continuation-ratio logistic model. This model compares each level of the ordinal variable to all lower levels. This can be modeled using binary logistic regression techniques, but new variables need to be constructed from the data set to allow this. The R package rms has a function cr.setup() which is a utility for preparing an outcome variable for a continuation ratio model.

I still fit the model to individual country to see where the issues occur. When the absolute Value of the predictor is greater than 10, I make table for that independent variable and the outcome for that country. ## Calculate the confident interval I use Wald test to obtain the CI because the likelihood ratio test doesn’t work

countryList <- c("Bangladesh", "Cameroon", "India", "Indonesia", "Kenya", "Nepal", "Philippines", "Uganda", "Zimbabwe")

fun <- function(data, countryList) {
  allData <- list()

  for (ctry in countryList) {
    m <- polr(formula = papm ~ gender + age + education + marital + occupation + live_area + mobile + internet, data = dat1 %>% filter(country == ctry), Hess = TRUE)

    ci <- confint.default(m)

    tmp <- as.data.frame(round(exp(cbind(OR = coef(m), ci)), 2)) %>%
      mutate(
        country = ctry,
        predictor = row.names(m)
      )

    allData[[ctry]] <- tmp
  }

  combinedData <- do.call(rbind, allData)

  return(combinedData)
}


tmp <- fun(dat1, countryList)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in polr(formula = papm ~ gender + age + education + marital + occupation
## + : design appears to be rank-deficient, so dropping some coefs
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in polr(formula = papm ~ gender + age + education + marital + occupation
## + : design appears to be rank-deficient, so dropping some coefs

## Warning in polr(formula = papm ~ gender + age + education + marital + occupation
## + : design appears to be rank-deficient, so dropping some coefs

Visualize the result

Clean the result

tmp1 <- tmp %>%
  mutate(independant_var = row.names(tmp)) %>%
  separate(independant_var, c("count", "var"), sep = "\\.") %>%
  mutate(var = str_replace_all(var, c(
    "gender" = "gender.",
    "age" = "age.Age",
    "education" = "education.",
    "marital" = "marital.",
    "occupation" = "occupation.",
    "live_area" = "live_area.",
    "mobile" = "mobile.",
    "internet" = "internet."
  ))) %>%
  separate(var, c("var", "level"), "\\.") %>%
  mutate(
    OR = ifelse(OR >= 1000, 0, OR),
    `2.5 %` = ifelse(`2.5 %` >= 1000, 0, `2.5 %`),
    `97.5 %` = ifelse(`97.5 %` >= 1000, 0, `97.5 %`)
  ) %>%
  mutate(level = as.factor(level))

Gender

tmp2 <- tmp1 %>%
  filter(var == "gender")

p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
  geom_line(group = 1, position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
  geom_point(size = 3, position = position_dodge(0.3)) +
  scale_y_continuous(trans = "log10") +
  scale_fill_brewer(palette = "Set1")+
  labs(title = "Reference - Male"
  )


ggplotly(p)

Age

tmp2 <- tmp1 %>%
  filter(var == "age") %>%
  mutate(level = as.factor(level))

p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
  geom_line(group = 1, position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
  geom_point(size = 3, position = position_dodge(0.3)) +
  scale_y_continuous(trans = "log10") +
  scale_fill_brewer(palette = "Set1")+
  labs(title = "Reference - Contunious"
  )


ggplotly(p)

Education

tmp2 <- tmp1 %>%
  filter(var == "education") %>%
  mutate(level = as.factor(level))

p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
  geom_line(group = 1, position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
  geom_point(size = 3, position = position_dodge(0.3)) +
  scale_y_continuous(trans = "log10") +
  scale_fill_brewer(palette = "Set1")+
  labs(title = "Reference - No formal education"
  )


ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Warning: `position_dodge()` requires non-overlapping x intervals

Marital

tmp2 <- tmp1 %>%
  filter(var == "marital") %>%
  mutate(level = as.factor(level))

p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
  geom_line(group = 1, position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
  geom_point(size = 3, position = position_dodge(0.3)) +
  scale_y_continuous(trans = "log10") +
  scale_fill_brewer(palette = "Set1")+
  labs(title = "Reference - Never married"
  )


ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Warning: `position_dodge()` requires non-overlapping x intervals

Occupation

tmp2 <- tmp1 %>%
  filter(var == "occupation") %>%
  mutate(level = as.factor(level))

p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
  geom_line(group = 1, position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
  geom_point(size = 3, position = position_dodge(0.3)) +
  scale_y_continuous(trans = "log10") +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "Reference - Employed"
  )


ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Warning: `position_dodge()` requires non-overlapping x intervals

Live area

tmp2 <- tmp1 %>%
  filter(var == "live_area") %>%
  mutate(level = as.factor(level))

p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
  geom_line(group = 1, position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
  geom_point(size = 3, position = position_dodge(0.3)) +
  scale_y_continuous(trans = "log10") +
  scale_fill_brewer(palette = "Set1")+
  labs(title = "Reference - Urban"
  )


ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Warning: `position_dodge()` requires non-overlapping x intervals

Mobile

tmp2 <- tmp1 %>%
  filter(var == "mobile") %>%
  mutate(level = as.factor(level))

p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
  geom_line(group = 1, position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
  geom_point(size = 3, position = position_dodge(0.3)) +
  scale_y_continuous(trans = "log10") +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "Reference - No"
  )


ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Warning: `position_dodge()` requires non-overlapping x intervals

Internet

tmp2 <- tmp1 %>%
  filter(var == "internet") %>%
  mutate(level = as.factor(level))

p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
  geom_line(group = 1, position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
  geom_point(size = 3, position = position_dodge(0.3)) +
  scale_y_continuous(trans = "log10") +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "Reference - No"
  )


ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

We can use anova() function to obtain the P value for each independent variable using likeli-hood ratio test